/* ========================================================================
 * 推文：资本存量计算那些事儿 (附各省、城市资本存量数据)
 * Data：CSMAR
 * Date：2022--2-16
 * Notes：以下代码的运行环境为 MacOS Stata17 MP，
 *        若在 Windows 环境下，应注意文件生成时的路径，
 *        并将斜杠 / 改为反斜杠 \，
 *        有任何问题，请邮件联系：chenbo2019@email.szu.edu.cn
 * ======================================================================= */


* # 1. 数据处理
********************
* ## 1.1 新建数据文件
********************

efolder, cd("~/Desktop/ChinaK")
efolder, cd("~/Desktop/ChinaK/Prov") ///
	sub(raw_data result_data picture)

global raw_data "~/Desktop/ChinaK/Prov/raw_data"
global result_data "~/Desktop/ChinaK/Prov/result_data"
global picture "~/Desktop/ChinaK/Prov/picture"

* ## 1.2 清洗 CSMAR 数据
************************

cd "$raw_data"
*- 固定资本形成总额
copy "https://gitee.com/arlionn/ChinaCapitalStock/raw/master/data/CRE_Gdp03.xlsx" ///
	 "CRE_Gdp03.xlsx", replace

import excel "CRE_Gdp03.xlsx", firstrow clear
labone, nrow(1 2)
drop in 1/2
rename (_all)(year provcode province invest)
destring year provcode invest, replace force
drop if provcode == 142
keep if inrange(year, 2000, 2017)
xtset provcode year
save "$result_data/invest_2000_2017.dta", replace

*- 固定资产投资价格指数
copy "https://gitee.com/arlionn/ChinaCapitalStock/raw/master/data/CRE_Pi04.xlsx" ///
	 "CRE_Pi04.xlsx", replace

import excel "CRE_Pi04.xlsx", firstrow clear
labone, nrow(1 2)
drop in 1/2
rename (_all)(year provcode province index)
destring year provcode index, replace force
keep if inrange(year, 2000, 2017)

xtset provcode year
tsfill, full
bys provcode: fillmissing province
xtset provcode year
save "$result_data/index_2000_2017.dta", replace

* ## 1.3 清洗张军数据
***********************

copy "https://gitee.com/arlionn/ChinaCapitalStock/raw/master/data/资本存量换算—张军.xlsx" ///
	 "资本存量换算—张军.xlsx", replace

import excel "资本存量换算—张军.xlsx", sheet("Sheet2") clear
nrow
keep province provcode _2000b
rename (_2000b)(zj_2000)
destring provcode zj_2000, replace force
drop if provcode == .
replace provcode = provcode * 10000
label var zj_2000 "2000年不变价资本存量 (张军)"
save "$result_data/stock_base_zj.dta", replace

* ## 1.4 清洗单豪杰数据
***********************

copy "https://gitee.com/arlionn/ChinaCapitalStock/raw/master/data/资本存量换算-单豪杰.xlsx" ///
	 "资本存量换算-单豪杰.xlsx", replace

import excel "资本存量换算-单豪杰.xlsx", firstrow clear
gen shj_2000 = invest_2000 * c2000_1952
label var shj_2000   "2000年不变价资本存量 (单豪杰)"
drop invest_2000 c2000_1952
save "$result_data/stock_base_shj.dta", replace

* ## 1.5 合并数据
***********************

use "$result_data/invest_2000_2017.dta", clear
merge 1:1 provcode year using "$result_data/index_2000_2017.dta", keep(1 3) nogen
merge m:1 provcode using "$result_data/stock_base_zj.dta", ///
	keepusing(zj) keep(1 3) nogen
merge m:1 provcode using "$result_data/stock_base_shj.dta", ///
	keepusing(shj) keep(1 3) nogen

* ## 1.6 计算西藏价格指数
***********************

	// 西藏的价格指数使用新疆和青海的均值 (单豪杰，2008)
inlist2 provcode, values(650000, 630000)
bys year inlist2: egen index2 = mean(index)
replace index2 = . if inlist2 == .
bys year: fillmissing index2
replace index = index2 if provcode == 540000
drop index2 inlist2

* ## 1.7 计算实际资本形成总额
***************************

xtset provcode year
gen deflator = 1 if year == 2000
bys provcode: replace deflator = index[_n] * deflator[_n-1] / 100 if mi(deflator)
gen real_invest = invest / deflator
save "$result_data/calculate.dta", replace

* # 2. 资本存量计算 (简易法)
**************************
* ## 2.1 张军法
******************** 
use "$result_data/calculate.dta", clear

*- 基期资本存量 (张军等，2000)
gen stock1 = real_invest / 0.1 if year == 2000
replace stock1 = real_invest + stock1[_n-1] * (1- 0.096) if mi(stock1)

* ## 2.2 单豪杰法
******************** 

*- 几何平均增长率
bys year: egen total_invest = sum(real_invest)
xtset provcode year
gen rate = d.total_invest / total_invest[_n-1]
egen growth = mean(rate) if inrange(year, 2001, 2005)
fillmissing growth
drop total_invest rate

gen stock2 = real_invest[_n+1] / (growth + 0.1096) if year == 2000
replace stock2 = real_invest + stock2[_n-1] * (1- 0.1096) if mi(stock2)

label var deflator    "平减指数"
label var real_invest "实际固定资本形成总额"
label var stock1      "资本存量 (张军法, K_1=I_1/0.1)"
label var growth      "投资增长率"
label var stock2      "资本存量 (单豪杰法, K_1=I_2/(g+delta))"
save "$result_data/capital_stock_1.dta", replace

*- 对比
collapse (sum)stock1 (sum)stock2, by(year)
for var stock*: replace X = X / 10000
gen gap = stock1 - stock2

tw 	(scatter stock1 year, c(l)) ///
	(scatter stock2 year, c(l) m(T)) ///
	(bar gap year, barw(0.7) yaxis(2)), ///
	xlabel(2000(2)2017) xtitle("") ytitle("") ///
	ylabel(-100(50)200 -50 " " -100 " ", nogrid) ///
	ylabel(0(20)60, axis(2)) ///
	legend(ring(0) rows(3) position(11)  ///
		order(1 "张军方法" 2 "单豪杰方法" 3 "资本存量差")) ///
	graphregion(fcolor(white))
graph export "$picture/prov_easy_way.png", width(1000) replace

* # 3. 资本存量计算 (复杂法)
**************************
* ## 3.1 合并四川、重庆
******************** 
use "$result_data/calculate.dta", clear

inlist2 provcode, values(500000, 510000)
bys year inlist2: egen invest2 = sum(real_invest)
replace invest2 = . if inlist2 == .
bys year: fillmissing invest2
replace real_invest = invest2 if provcode == 510000
drop if provcode == 500000
drop invest2 inlist2
xtset provcode year

* ## 3.2 张军法
******************** 
gen stock3 = zj_2000 if year == 2000
replace stock3 = real_invest + stock3[_n-1] * (1- 0.096) if mi(stock3)

* ## 3.3 单豪杰法
******************** 
gen stock4 = shj_2000 if year == 2000
replace stock4 = real_invest + stock4[_n-1] * (1- 0.1096) if mi(stock4)

label var deflator    "平减指数"
label var real_invest "实际固定资本形成总额"
label var stock3      "资本存量 (张军法, K_1为1952年)"
label var stock4      "资本存量 (单豪杰法, K_1为1952年)"
save "$result_data/capital_stock_2.dta", replace

use "$result_data/capital_stock_1.dta", clear
merge 1:1 provcode year using "$result_data/capital_stock_2", keep(1 3) nogen
save "$result_data/capital_stock_prov.dta", replace
erase "$result_data/capital_stock_1.dta"
erase "$result_data/capital_stock_2.dta"

collapse (sum)stock1 (sum)stock2 (sum)stock3 (sum)stock4, by(year)
for var stock*: replace X = X / 10000
gen gap = stock3 - stock4

	// 1952年计价的比较
tw 	(scatter stock3 year, c(l)) ///
	(scatter stock4 year, c(l) m(T)) ///
	(bar gap year, barw(0.7) yaxis(2)), ///
	xlabel(2000(2)2017) xtitle("") ytitle("", axis(2)) ///
	ylabel(, nogrid) ylabel(0(5)20, axis(2)) ///
	legend(ring(0) rows(3) position(11) ///
		order(1 "张军方法 Base1952" 2 "单豪杰方法 Base1952" 3 "资本存量差")) ///
	graphregion(fcolor(white))
graph export "$picture/prov_ase_1952.png", width(1000) replace

local color1 = `"lcolor("240 59 32") mcolor("240 59 32")"'
local color2 = `"lcolor("49 163 84") mcolor("49 163 84")"'
local color3 = `"lcolor("254 178 76") mcolor("254 178 76")"'
local color4 = `"lcolor("161 217 155") mcolor("161 217 155")"'

	// 四种方法比较
tw 	(scatter stock1 year, c(l) m(oh) `color1') ///
	(scatter stock2 year, c(l) m(Th) `color2' lp(dash)) ///
	(scatter stock3 year, c(l) m(Sh) `color3') ///
	(scatter stock4 year, c(l) m(Dh) `color4' lp(dash)), ///
	xlabel(2000(2)2017) xtitle("") ///
	ylabel(, nogrid angle(0)) ///
	legend(ring(0) rows(4) position(11) ///
		order(1 "张军方法" 2 "单豪杰方法" 3 "张军方法 Base1952" 4 "单豪杰方法 Base1952")) ///
	graphregion(fcolor(white)) plotregion(margin(left))
graph export "$picture/prov_four_way.png", width(1000) replace



* # 4. 城市资本存量计算
**************************
* ## 4.1 新建数据文件
********************

efolder, cd("~/Desktop/ChinaK/City") ///
	sub(raw_data result_data picture)

global raw_data "~/Desktop/ChinaK/City/raw_data"
global result_data "~/Desktop/ChinaK/City/result_data"
global picture "~/Desktop/ChinaK/City/picture"


* ## 4.2 城市固定资产投资数据
*******************************
cd "$raw_data"
copy "https://gitee.com/arlionn/ChinaCapitalStock/raw/master/data/CRE_Invcsct.xlsx" ///
	 "CRE_Invcsct.xlsx", replace

import excel "CRE_Invcsct.xlsx", firstrow clear
labone, nrow(1 2)
drop in 1/2
drop Cttyp
drop if Ctnm == "巢湖市"
rename (_all)(year city citycode provcode province invest)
destring year *code invest, replace force
keep if inrange(year, 2000, 2017)


*- 整合平衡面板
bys citycode: drop if _N <= 10
xtset citycode year
tsfill, full
bys citycode: fillmissing city
bys citycode: fillmissing province
bys citycode: fillmissing provcode
xtset citycode year
order city citycode year province provcode
for var _all: cap format %15s X

*- 缺失值处理
replace invest = 112746900 if citycode == 120000 & year == 2017
replace invest = 1322000   if citycode == 131100 & year == 2001
replace invest = 14908000  if citycode == 150100 & year == 2017
replace invest = 29587900  if citycode == 150200 & year == 2017
replace invest = 1789100   if citycode == 150300 & year == 2017
replace invest = 15069500  if citycode == 150400 & year == 2017
replace invest = 14652100  if citycode == 150500 & year == 2017
replace invest = 30657415  if citycode == 150600 & year == 2017
replace invest = 391337    if citycode == 150700 & year == 2000
replace invest = 7090600   if citycode == 150700 & year == 2017
replace invest = 289501    if citycode == 150800 & year == 2000
replace invest = 328854    if citycode == 150800 & year == 2001
replace invest = 383500    if citycode == 150800 & year == 2002
replace invest = 6502100   if citycode == 150800 & year == 2017
replace invest = 226684    if citycode == 150900 & year == 2000
replace invest = 193583    if citycode == 150900 & year == 2001
replace invest = 315700    if citycode == 150900 & year == 2002
replace invest = 5334700   if citycode == 150900 & year == 2017
replace invest = 922951    if citycode == 450800 & year == 2004
replace invest = 294600    if citycode == 451000 & year == 2000
replace invest = 71504     if citycode == 451100 & year == 2000
replace invest = 63208     if citycode == 451100 & year == 2001
replace invest = 405500    if citycode == 451200 & year == 2000
replace invest = 180100    if citycode == 451300 & year == 2000
replace invest = 80099     if citycode == 451400 & year == 2000
replace invest = 111338    if citycode == 451400 & year == 2001
replace invest = 38506000  if citycode == 520100 & year == 2017
replace invest = 16525000  if citycode == 520200 & year == 2017
replace invest = 25236200  if citycode == 520300 & year == 2017
replace invest = 7945400   if citycode == 520400 & year == 2017
replace invest = 214411    if citycode == 530600 & year == 2000
replace invest = 133224    if citycode == 530700 & year == 2001
replace invest = 125938    if citycode == 530800 & year == 2000
replace invest = 213654    if citycode == 530800 & year == 2001
replace invest = 242230    if citycode == 530900 & year == 2000
replace invest = 201661    if citycode == 530900 & year == 2001
replace invest = 25898800  if citycode == 610300 & year == 2015
replace invest = 200900    if citycode == 611000 & year == 2000
replace invest = 151765    if citycode == 620600 & year == 2000
replace invest = 183858    if citycode == 620700 & year == 2000
replace invest = 233794    if citycode == 620700 & year == 2001
replace invest = 298873    if citycode == 620800 & year == 2000
replace invest = 213968    if citycode == 620800 & year == 2001
replace invest = 239007    if citycode == 620900 & year == 2000
replace invest = 350012    if citycode == 620900 & year == 2001
replace invest = 204675    if citycode == 621000 & year == 2000
replace invest = 257512    if citycode == 621000 & year == 2001
replace invest = 152718    if citycode == 621100 & year == 2000
replace invest = 167587    if citycode == 621100 & year == 2001
replace invest = 209774    if citycode == 621100 & year == 2002
replace invest = 133458    if citycode == 621200 & year == 2000
replace invest = 143789    if citycode == 621200 & year == 2001
replace invest = 174038    if citycode == 621200 & year == 2002
replace invest = 221297    if citycode == 621200 & year == 2003
replace invest = 114000    if citycode == 640400 & year == 2001

xtset citycode year
by citycode: mipolate invest year, nearest gen(invest_mi)
replace invest = invest_mi if mi(invest)
drop invest_mi

label var invest "固定资产投资 (万元)"
save "$result_data/invest_city.dta", replace

* ## 4.3 固定资产投资价格指数
*****************************
copy "https://gitee.com/arlionn/ChinaCapitalStock/raw/master/data/CRE_Pi04.xlsx" ///
	 "CRE_Pi04.xlsx", replace

import excel "CRE_Pi04.xlsx", firstrow clear
labone, nrow(1 2)
drop in 1/2
rename (_all)(year provcode province index)
destring year provcode index, replace force
keep if inrange(year, 2000, 2017)

xtset provcode year
tsfill, full
bys provcode: fillmissing province
xtset provcode year
save "$result_data/index_2000_2017.dta", replace

* ## 4.4 计算资本存量
*****************************

*- 合并数据
use "$result_data/invest_city.dta", clear
merge m:1 provcode year using "$result_data/index_2000_2017.dta", ///
	keep(1 3) nogen

*- 计算实际投资
xtset citycode year
gen deflator = 1 if year == 2000
bys citycode: replace deflator = index[_n] * deflator[_n-1] / 100 if mi(deflator)
gen real_invest = invest / deflator

*- 张军法
gen stock1 = real_invest / 0.1 if year == 2000
replace stock1 = real_invest + stock1[_n-1] * (1- 0.096) if mi(stock1)

*- 单豪杰法
bys year: egen total_invest = sum(real_invest)
xtset citycode year
gen rate = d.total_invest / total_invest[_n-1]
egen growth = mean(rate) if inrange(year, 2001, 2005)
fillmissing growth
drop total_invest rate

gen stock2 = real_invest[_n+1] / (growth + 0.1096) if year == 2000
replace stock2 = real_invest + stock2[_n-1] * (1- 0.1096) if mi(stock2)
drop growth

label var deflator    "平减指数 (2000=1)"
label var real_invest "实际固定投资 (2000不变价)"
label var stock1      "资本存量 (张军法, 单位：万元)"
label var stock2      "资本存量 (单豪杰法, 单位：万元)"
save "$result_data/capital_stock_city.dta", replace

use "$result_data/capital_stock.dta", clear
collapse (sum)stock1 (sum)stock2, by(year)
for var stock*: replace X = X / 100000000
gen gap = stock1 - stock2

tw 	(scatter stock1 year, c(l)) ///
	(scatter stock2 year, c(l) m(T)) ///
	(scatter gap year, m(oh) mcolor(blue) yaxis(2)), ///
	xlabel(2000(2)2017) xtitle("") ///
	ytitle("资本存量 (万亿)") ytitle("资本存量差 (万亿)", axis(2)) ///
	ylabel(0(40)220, nogrid) ///
	ylabel(6(2)14, axis(2)) ///
	legend(ring(0) rows(3) position(5)   ///
		order(1 "张军方法" 2 "单豪杰方法" 3 "资本存量差")) ///
	graphregion(fcolor(white)) 
graph export "$picture/city_easy_way.png", width(1000) replace
